arXiv:1509.03168vl [astro-ph.HE] 10 Sep 2015 


Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 11 September 2015 (MN OTeX style file v2.2) 


Three-dimensional simulations of super-critical black hole accretion 
disks — luminosities, photon trapping and variability. 

Aleksander S^dowski^* and Ramesh Narayan^* 

^ MIT Kavli Institute for Astrophysics and Space Research 77 Massachusetts Ave, Cambridge, MA 02139, USA 
^ Harvard-Smithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA 02134, USA 


11 September 2015 


ABSTRACT 

We present a set of four three-dimensional, general relativistic, radiation MHD simula¬ 
tions of black hole accretion at super-critical mass accretion rates, M > MEdd- We use these 
simulations to study how disk properties are modified when we vary the black hole mass, 
the black hole spin, or the mass accretion rate. In the case of a non-rotating black hole, we 
find that the total efficiency is of order 3%Mc^, approximately a factor of two less than the 
efficiency of a standard thin accretion disk. The radiation fiux in the funnel along the axis is 
highly super-Eddington, but only a small fraction of the energy released by accretion escapes 
in this region. The bulk of the 3%Mc^ of energy emerges farther out in the disk, either in the 
form of photospheric emission or as a wind. In the case of a black hole with a spin parameter 
of 0.7, we find a larger efficiency of about 8%Mc^. By comparing the relative importance 
of advective and diffusive radiation transport, we show that photon trapping is effective near 
the equatorial plane. However, near the disk surface, vertical transport of radiation by diffu¬ 
sion dominates. We compare the properties of our fiducial three-dimensional run with those 
of an equivalent two-dimensional axisymmetric model with a mean-field dynamo. The latter 
simulation runs nearly 100 times faster than the three-dimensional simulation, and gives very 
similar results for time-averaged properties of the accretion flow. 

Key words: accretion, accretion discs - black hole physics - relativistic processes - methods: 
numerical 


1 INTRODUCTION 

Black hole (BH) accretion disks are common in the Universe. It 
appears that virtually every galaxy harbours a supermassive black 
hole (SMBH) in its nucleus and it is likely that every one of these 
SMBHs has some kind of an accretion flow. Moreover, just in our 
own Galaxy, there are probably thousands of stellar-mass BHs in 
binaries accreting gas from their companions, of which a few dozen 
have been detected in X-rays and are widely studied. 

Because of the compactness of BHs, accreting gas can liber¬ 
ate significant amounts of gravitational energy and make the ac¬ 
creting systems extraordinary luminous. Moreover, magnetic fields 
near the BH encourage the extraction of rotational energy of spin¬ 
ning BHs, leading to formation of powerful relativistic jets. 

Early theoretical work on accretion disks was limited to one¬ 
dimensional analytical models. Later, 1 + 1 and two-dimensional 
models with cr-viscosity were developed. Because accretion flows 
are by nature turbulent, such simplified models were not adequate 
in most cases. Moreover, they often made strong assumptions, e.g., 
no mass loss, constant a, zero-torque at the innermost stable circu- 
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lar orbit (ISCO), which may not be satisfied in real systems. This 
motivated the development of techniques for numerically modeling 
multi-dimensional turbulent accretion flows. 


The key breakthroughs were the identification of the mag- 
netorotational instability (MRI, [Balbus & Hawley 1 199 l| l and the 
development of magnetohydrodynamical (MHD) numerical codes, 
both Newtonian and relativistic. Initial efforts were focused on sim¬ 
ulating optically thin disks, as are likely to be present at the low¬ 
est accretion rates. Such systems are relatively simple to simulate 
because the radiation is weak and is, moreover, decoupled from 
the gas. Only in recent years have more advanced radiation-MHD 
(RMHD) codes been developed which can be used to study radia- 
tively luminous systems. The pioneering initial work was based on 
Newtonian codes with crude (flux-limited diffusion) radiative trans¬ 
port { [Ohsuga et al.||200^ |Ohsuga & Mineshi^|2011| l. This was 
later followed by codes using more advanced radiative transport 
schemes, either still in the Newtonian or special-relativistic approx¬ 
imation ( [Jiang et ^ 2012||2014a[ Ohsuga & Takahashi|2015jl, < 


in general relativity (Sadowski et al.|20131|McKinney et al.|20141 
[Fragile et al.|2014t . 

So far global simulations of optically thick disks that evolve 
the radiation field in parallel to gas have been performed only for 
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super-critical (exceeding the Eddington value, see equation ac¬ 
cretion rates Such disks are geometrically thick and do not re¬ 
quire excessive resolution near the equatorial plane which, so far, 
makes self-consistent simulations of thinner disks too expensive. 

In this paper we continue the numerical study of super-critical 
BH accretion flows by performing a set of four three-dimensional 
simulations, the parameters of which probe different accretion 
rates, BH spins, and BH masses. Such a comparative study us¬ 
ing a single code (and set of assumptions) has not yet been per¬ 
formed. In addition, we compare the properties of our fiducial 
three-dimensional model with an equivalent axisymmetrical two- 
dimensional run which is simulated using the mean-field dynamo 
model described in |S^dowski et al.| ( [2015a| ). 

We begin with a description of the numerical methods (Sec¬ 
tion and details of the five simulations (Section [^. We then dis¬ 
cuss the results, focusing on the luminosities of the simulated disks 
(Section 1^, the efficiency of photon trapping (Section]^, and the 
variability of the emitted radiation (Section [^. Finally, we assess 
the strengths and weaknesses of 2D simulations (Section [^, and 
list the conclusions in the Summary (Section [^. 


2 NUMERICAL METHODS 


The simulations described in this paper were performed with the 
general relativistic radiation magnetohydrodynamical (GRRMHD) 
code KORAL ( [S^dowski et al.||2013| ) which solves the conserva¬ 
tion equations in a fixed, arbitrary spacetime using finite-difference 
methods. The equations we solve are. 


(pm");;, = 0, 

(1) 

II 

g 

(2) 

= -G„ 

(3) 


where p is the gas density in the comoving fluid frame, are the 
components of the gas four-velocity as measured in the “lab frame”, 
Ty is the MHD stress-energy tensor in this frame, 

= (p + Ug + pg + b^)u^Uy + (pg + “ b^by, (4) 

Ry is the stress-energy tensor of radiation, and Gy is the radiative 
four-force describing the interaction between gas and radiation (see 


S^dowski et al. 2014 for a more detailed description). Here, Ug and 


Pg = (T - l)Ug represent the internal energy and pressure of the 
gas in the comoving frame and b^ is the magnetic field 4-vector 
( [Gammie et al.|[200^ . The magnetic pressure is p^ag = in 
geometrical units. 

The magnetic field is evolved via the induction equation. 


dti ^PgB') = -dj ( - b'u^)) , 


(5) 


where is the magnetic field three-vector ( [k 


.omissarov 


1999| ), and 


is the metric determinant. The divergence-free criterion is en¬ 
forced using the flux-constrained scheme of |T6th| ( [2000| ). 

The radiation field is evolved through its energy density and 
flux, and the radiation stress-energy tensor is closed by means of the 
Ml closure scheme ( |Levermore|1984[rS^dowski et al.|201^ . The 
energy exchange between gas and radiation is by free-free emis¬ 
sion/absorption as well as Compton scattering. The latter is treated 


^ A number of groups (e.g., |Shafee et al.|2008a|[^hnittman et al.|2Q13[ 
[Avara et al.|2015) have performed simulations of thin disks in pure hydro- 
dynamical setup implementing arbitrary cooling function. 


in the “blackbody” Comptonization approximation as described in 
[S^dowski & Naray^ ( |2015c] ). 

Four of the five simulations described here were performed in 
three dimensions (3D), while the fifth was carried out in 2D, as¬ 
suming axisymmetry and using the mean-field dynamo model de¬ 
scribed in |S^dowski et al.|(2015a| ) with model parameters identical 
to those used there. 

We use modified Kerr-Shild coordinates with the inner edge of 
the domain inside the BH horizon. The simulations were run with a 
moderately high resolution of 252 grid cells spaced logarithmically 
in radius, 234 grid cells in the polar angle, concentrated towards the 
equatorial plane, and 128 cells in azimuth. 

All details of the numerical method are given in|S^dowski et| 

|^{ [20T4l ). 


3 SIMULATIONS 


3.1 Units 


We adopt the following definition for the Eddington mass accretion 
rate. 


- 



( 6 ) 


where L^dd = 1-25 x lO^^M/M© ergs/s is the Eddington luminosity, 
7] is the radiative efficiency of a thin disk around a black hole with 
a given spin = a/M, 


ri=\- 



(7) 


and risco is the radius of the Innermost Stable Circular Orbit 
(ISCO). According to this definition, a standard thin, radiatively 
efficient dis k ([Shakura & Sunyaev| 197^ [Novikov & Thorne|1973[ 

Frank et ^|1985| ) accreting at Mndd would have a luminosity of 

L^dd- For zero BH spin, MEdd = 2.48 X lO^^M/M© g/s. 

Hereafter, we use the gravitational radius Vg = GM/d as the 
unit of length, and r^/c as the unit of time. 


3.2 Initial setup 

Each of the five simulations was initialized as in ISadowski et al.l 
P014| t, viz., by starting with the hydrodynamical equilibrium torus 
of |Penna et alJh[2013a| ) with the angular momentum parameters 
listed in Table and then redistributing the pressure between gas 
and radiation such that local thermal equilibrium is established with 
the total pressure unchanged. The resulting, radiation-pressure sup¬ 
ported torus is close to equilibrium. The initial density was set 
through the torus entropy parameter and was adjusted to provide an 
optically thick torus that would give super-critical accretion once 
the simulation has reached steady state. 

The initial torii were threaded by weak poloidal magnetic 
field. As each simulation proceeded, the field grew in strength and 
led to the onset of the magnetorotational instability, which triggered 
and maintained MHD turbulence in the disk. For most models we 
started with multiple loops with alternating polarity. For one model 
(rOll) we used a single loop. Both field configurations were ini¬ 
tialized as in |S^dowski et al.| ( [2015a| ). 

We performed two simulations with BH mass lOM© and zero 
BH spin. One of these (run rOOl) resulted in a mean accretion rate 
of ~ lOMEdd. while the other (r003), which was initialized with 
a more optically thick torus, accreted at ~ llSM^dd- These two 
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models allowed us to study the role of the mass accretion rate on 
disk properties. 

The third model (rOll) was initialized with a torus similar 
to model rOOl, but we assumed a rotating BH with spin param¬ 
eter a* = 0.7. This model was the only one of the five that was 
initialized with a single poloidal loop of magnetic field. The hope 
was that the single loop would lead to a strong magnetic field at 


the horizon and would give a magnetically arrested disk ( Igumen 


shchev et al.|2003[|Narayan et al.|2003[|Tchekhovskoy et al.|2011[ 
McKinney et al.|2012| l. However, although this simulation was run 
up to a time of nearly 15,000, this duration was still insufficient 
to reach the MAD limit. Therefore, the magnetic field at the BH 
still did not reach the saturated level appropriate to the MAD state. 
Nevertheless, by comparing this model with rOOl, we were able to 
study the effect of BH spirj^ 

In the fourth model (r020), we increased the BH mass to 
IOOOMq BH, but kept the mass accretion rate at ~ lOMEdd- This 
enabled us to investigate the role of BH mass. 

All of the above models were run in 3D. The fifth model 
(d300) was a 2D axisymmetric simulation which used the mean- 
field dynamo model of |S^dowski et al.| ( |2015a| ). This model was 
initialized with exactly the same torus configuration as in model 
rOOl. The only difference was that it was evolved in 2D instead 
of 3D. The purpose of this model was to assess how well the 2D 
dynamo model captures the properties of the 3D model. 


3.3 Accretion flow properties 

Each of the four 3D simulations was run up to a final time tmax ~ 
15,000 - 20,000, by which time all had developed quasi-steady, 
turbulent accretion via optically and geometrically thick disks. The 
time histories of the mass accretion rate through the BH horizon are 
shown in Fig.[2 In all the runs, gas starts crossing the BH horizon 
at a substantial rate only after t ^ 2000. This is the time needed for 
the magnetorotational instability to make the disk turbulent, and 
for gas from the inner edge of the initial torus to accrete on the BH. 
Once accretion begins, the mass accretion rate increases rapidly. 
In fact, M overshoots and remains somewhat enhanced until time 
t ^ 9000, and only beyond this time does the accretion rate settle 
down to a quasi-steady value. In the following, we focus on disk 
properties averaged over the latter quasi-steady stage of accretion, 
from t - 9000 iot - tmax- 

The radial profile of the mass accretion rate is given by 

pn p2n 

M = - I V^pw"d0d6>. (8) 

Jo Jo 

Fig. shows the time-averaged profiles of this quantity as a func¬ 
tion of radius r corresponding to the four 3D runs. The flat sections 
of the profiles at relatively small radii denote the region where the 
flow has reached inflow equilibrium. Given the somewhat limited 
duration of the simulations, and the relatively low radial velocity 
of the accreting gas in the disk, inflow equilibrium is reached only 
up to a radius r^q ~ 20 - 30. The outflows in the jet and wind re¬ 
gions have larger velocities, and therefore the outflowing regions 
are causally connected with the equatorial disk out to much larger 
distances from the BH. In particular, the funnel region, which is 
filled in most cases with gas escaping at u > 0.1c, achieves equilib¬ 
rium all the way out to the domain boundary at rout = 1000. 

Figure shows poloidal and azimuthal slices through model 


In a recent paper, [^Kinney et~n^ ( |2015t present true MAD simulations 


0 0.2465 0.493 0.7395 0.986 
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Figure 1. Time history of the mass accretion rate at the BH in Eddington 
units (see §3.1 [ for definition) for the four three-dimensional models consid¬ 
ered in this paper. Model parameters are given in Table[^ 



Figure 2. Time- and azimuth-averaged radial profiles of the mass accretion 
rates (eq.jSj in the four three-dimensional simulations. 


rOOl at time t = 16400. The colors in the left halves of the panels 
show the magnitude of the radiation field, while those in the right 
halves show the gas density. The arrows in the two halves show the 
direction of the radiative flux and gas velocity, respectively, with 
the arrow thicknesses indicating the corresponding magnitudes. For 
clarity, the vector fields were computed from time-averaged output. 

The gas is concentrated near the equatorial plane and forms a 
geometrically thick (or slim) disk with density scale-height HjR ^ 
0.25. Non-axisymmetric modes are clearly visible in the right 
panel, showing the value of running the simulations in 3D. Because 
of the large density, the gas is optically thick and advects with it 
a significant amount of radiation. This explains why the radiation 
flux has its largest magnitude in the bulk of the disk. 

The gas in the disk region around the mid-plane is turbulent. 
On average the gas moves there slowly towards the BH. Outside the 
bulk of the disk, within ~ 40 deg from the pole, the gas flows out¬ 
ward, being driven mostly by the radiation pressure force exerted 
by the radial radiative fiu}|^ 

Gas in the disk region orbits around the BH but because the 


in the radiative super-critical regime. Those represent a different class of 
accretion flows than the ones considered here. 

^ For a comprehensive study of the acceleration mechanism of outflows in 
hydrodynamical and radiative disks see |Moller & S^dowski| ( |2015j . 
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Table 1. Model parameters 



rOOl 

r003 

rOll 

r020 

d300 (2D) 

Mbu 

lOMo 

lOMo 

lOMo 

lOOOMo 

lOMo 

MIM^dd 

10.0 

175.8 

9.7 

10.1 

8.9 


0.0 

0.0 

0.7 

0.0 

0.0 

P0,init 

4.27 X 10-3 

6.61 X 10-3 

4.68 X 10-3 

3.99 X 10-3 

4.27 X 10-3 

/^max 

10.0 

10.0 

10.0 

10.0 

10.0 

initial B loops 

multi. 

multi. 

poloidal 

multi. 

multi. 

Nr X No X N^ 

252 X 234 X 128 

252x234x 128 

252 X 234 X 128 

252 X 234 X 128 

252 X 234 X 1 

r min / r max / ^0 / Hq 

1.85/ 1000/0/0.6 

1.85/ 1000/0/0.6 

1.85/ 1000/0/0.6 

1.85/1000/0/0.6 

1.85/ 1000/0/0.6 

^1 n! r 2 ln^ 

0.705/40/1000/22.5 

0.705/40/1000/22.5 

0.705/40/1000/10 

0.705/40/1000/22.5 

0.705/40/1000/22.5 

^max 

20,000 

19,000 

16,100 

19,200 

190,000 


M - average accretion rate 

Po,init - maximal density of the initial torus in g/cm^ 

/^max - maximal value of initial total to magnetic pressure ratio 
Nr X No X N(i) - resolution _ 


^min / rmax 1 fQ 1 Hq - grid parameters defined in 

Sadowski et al. J2015a 


^ / ri / r 2 / Tin - parameters of the initial torus as defined in 

Penna et al. 

j2013a| 


ftnay " durution of simulation in units of GMjc^ 



Figure 4. Average radial profiles of the density-weighted gas angular mo¬ 
mentum, Ufj). The black lines show the Keplerian profiles for spin = 0.0 
(solid) and a* = 0.7 (dashed). 


disk is geometrically thick the rotation is mildly sub-Keplerian, as 
shown by Fig. Outside the ISCO, the deviation from the Keple¬ 
rian profile is no more than 13% in any of the models. The angu¬ 
lar momentum profile flattens towards the BH horizon, reflecting 
the fact that the effective visosity, which transports angular mo¬ 
mentum, is less effective in the plunging region. The profile is flat¬ 
test for the models with the lowest accretion rate. Extrapolating to 
sub-Eddington accretion rates, we thus expect the viscous torques 


to become insignificant for thin disks, in agreement with previous 

work[Paczynski (2000J; 

Afshordi & Paczynski[([2003[);[Shafee et al.[ 

([2008 a|b[); [Penna et al.[( 

2oTor 


The top set of panels in Eig. shows the time- and azimuth- 
averaged distribution of density for all the five simulations. Stream¬ 
lines reflect the velocity held of the gas, with the thickness of lines 
denoting the density-weighted average velocity. Eor all the simula¬ 
tions the accretion flow is geometrically thick with density peaking 
at the equatorial plane. The density values are similar in all the 
runs except the simulation with the highest accretion rate, r003, 
which has a significantly higher density of gas (and increased op¬ 
tical depth). In all cases the gas flows on average towards the BH 
deep in the disk. The velocity of gas in the funnel is pointing out¬ 


ward outside the stagnation radius separating the polar region of in¬ 
flowing gas (which is the boundary condition imposed by the pres¬ 
ence of BH horizon) and the outflowing gas, driven either by radia¬ 
tion pressure or pressure gradients. In the case of simulations with 
non-rotating BHs and moderate accretion rates (three-dimensional 
models rOOl and r020), the stagnation radius at the axis is located 
near r ~ 10. Eor the run with much higher accretion rate (r003) it 
shifts significantly outward to r ~ 50 due to the much larger opac¬ 
ity, which prevents radiation from ejecting gas from the innermost 
region. On the other hand, the stagnation radius is very close to the 
BH horizon for the simulation with a rotating BH. In this model, 
gas is accelerated outward by an additional energy source, viz., the 
spin energy of the BH. The transition between the region of inflow 
(inside the disk), and outflow (in the funnel) is quite rapid — when 
gas particles are blown out of the disk and enter the funnel they 
quickly gain large outward radial velocity and join the outflow. 

The bottom set of panels in Eig. shows the averaged radia¬ 
tive flux for the five simulations. The colors denote the magnitude 
and arrows show the direction of the radiative flux. The radiation 
emitted by hot gas in the disk is mostly advected with the optically 
thick gas, but in regions close to the disk surface it also diffuses 
down the local density gradient ( [Jiang et ^|2014b| ). As a result, 
some fraction of the radiation is advected into the BH and the rest 
naturally Alls up the funnel region. There, radiation pressure ac¬ 
celerates gas outward along the axis of the funnel, resulting in a 
jet that travels at a modest fraction (0.2-0.5) of the speed of light 
( jS^dowski & Narayan|2015b] t. The properties of the radiation held 
depend on the accretion rate and BH spin. Eor model r6)03, which 
has the largest accretion rate, a significantly higher fraction of all 
the radiative flux in the inner region is advected into the BH, even 
in the inner part of the funnel region. In simulation rOl 1, the rota¬ 
tional energy of the BH is extracted (although not very efficiently 
due to the sub-MAD level of magnetic flux accumulated at the hori¬ 
zon) through the Blandford-Znajek process. This extra energy is 
converted into radiation and therefore this model has the highest 
radiative flux magnitude in the funnel region among all the runs. 

The qualitative properties of the 3D simulations described 
here agree well with accretion flows simulated in the recent years 
in axisymmetry e.g., by Sqdowski et al.j ( |2014j ), and in three- 
dimensions by [McKinney et al.| ( |2014j |. We comment on the differ¬ 
ences between our models and the simulation presented by [Jiang[ 
[et al.[ ( [2014b| ) in Section [STT] Below we discuss in detail three as- 
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poloidal plane 


equatorial plane 



Figure 3. Shows the magnitude of the radiative flux (orange colors in the left 
panel) for a snapshot taken near the end of the rOO 1 simulation. The left and 
respectively. Streamlines of the radiative flux and gas velocity are azimuth- and 
of the respective quantity. 



_50 -40 -30 -20 -10 0 10 20 30 40 50 


r/rg 

half of each panel) and the gas density (grey colors in the right half of each 
right panels correspond to slices through the poloidal and equatorial planes, 
time-averaged. The thickness of the streamlines increases with the magnitude 


pects of our models - luminosities, efficiency of photon trapping, 
and variability. 


4 LUMINOSITIES 

Accretion takes place if there is an efficient mechanism for trans¬ 
porting angular momentum outward. In BH accretion disks, we 
believe that this transport results from turbulence sustained by 
magnetorotational-unstable magnetic field. The angular momen¬ 
tum flux is accompanied by fluxes of energy in various forms. In the 
standard picture of a thin disk, the accreting gas brings in kinetic 
(orbital and turbulent), thermal, and binding energy. The latter has 
a negative sign, so effectively the flux of binding energy transports 
energy out (and deposits it at infinity). The exchange of angular 
momentum would not be possible without a shear stress (“viscos¬ 
ity”) which again causes an outward flux of energy. For turbulent 
magnetic disks, this energy flux comes from the work done by mag¬ 
netic fields. In addition, there is radiation flux. In a thin disk, radia¬ 
tion carries energy out to infinity. However, in geometrically thick 
super-Eddington accretion flows, the radiation energy flux can be 
either outward or inward, depending on how effectively the radia¬ 
tion is trapped in the optically thick flow. In addition, these models 
can also have mechanical energy flowing out in a wind or jet. 

We postpone a comprehensive discussion of the energetics in 
multi-dimensional accretion disk to a forthcoming paper. Below, 
we limit ourselves to simple angle-integrated luminosities of the 
simulated disks. 


4.1 Total luminosities 


The most fundamental definition of the luminosity is. 


■ r r(r[+R^+pu^)d=gdm 

Jo Jo 


(9) 


which is the total rate of energy flowing through a sphere at some 
radius r at some instant of time. This quantity accounts for all forms 
of energy except the rest-mass energy (which has been subtracted 



r/M 


Figure 5. Total luminosity (solid lines, eq. and radiative luminosity 
(dashed lines, eq. ED integrated over the whole sphere for the four three- 
dimensional simulations. 


out via the term pu^). In particular, Ltotai includes binding (gravi¬ 
tational), radiative, kinetic, thermal and magnetic energies. In an 
equilibrium steady state accretion disk, the the time-averaged lu¬ 
minosity Ltotai is conserved, i.e., it has the same value at all radii. 
This is because of energy conservation: any energy that appears in 
one form must ultimately come from one of the other forms dis¬ 
cussed above, and therefore the sum of all forms of energy remains 
constant. 

The radial profile of total luminosity is plotted with solid lines 
in Fig. We see that it is indeed constant to good accuracy for 
all the runs. All the simulations with non-rotating BHs have to¬ 
tal luminosity close to 3%Mc^. For an accretion rate of lOMEdd 
this corresponds to ~ 5LEdd (for the adopted definition of MEdd see 
equation®. The luminosity is significantly higher ~ 8%Mc^ for 
the simulation (rOll) with a rotating BH. There are two reasons 
for this: (i) the disk around a spinning BH extends deeper into the 
potential well since the ISCO is at a smaller radius (correspond- 
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ingly, the thin disk efficiency for a, = 0.7 is //thin ~ 10.3%Mc^), and 
(ii) the accumulated magnetic flux at the BH horizon allows for the 
extraction of BH kinetic energy through the Blandford-Znajek pro¬ 
cess. The predicted rate of that process (for the measured average 
magnetic flux parameter 0^15 and = 0.7) is y/jet ~ 6%. The to¬ 
tal energy available is therefore ^thin+^jet ~ 16%. In actuality, ~ 8% 
is extracted, which is roughly 50% of the total available. A similar 
50% energy extraction is seen also in the case of non-rotating BHs: 
a thin disk around a non-spinning BH has an efflciency of 5.7% 
whereas, as mentioned above, the luminosity of our models is only 
3%. Total luminosities for all the simulations are given in Edding¬ 
ton units in the third row of Table [2] 

The total luminosity as deflned above may be, in principle, 
sensitive to the initial conditions. In the ideal world, one would ini¬ 
tiate a simulation by dumping marginally bound (zero Bernoulli 
number) gas from inflnity. Because the duration of real simula¬ 
tions is limited, we have to start our simulatons from a bound torus 
located close to the BH with its inner edge at rin = 22.5. The 
Bernoulli number, deflned as 


Be = - 


Tj- + pu! 
pu^ 


( 10 ) 


of the initial gas at the very inner edge is Be ^ -0.014 and ap¬ 
proaches zero inversely proportional to radius. Thus, in principle, 
the luminosity estimates given above may be overestimated by 
up to ~ l%Mc^ (it cannot be more, but it may be less if there 
is signiflcant radial mixing of gas in the initial torus). To esti¬ 
mate how strongly the initial Be affects results we compared the 
total luminosities averaged over t = 12,000 - 13,000 and t = 
18,000 - 19,000. At t ~ 12,000, the BH had accreted an amount 
of gas equivalent to that contained in the initial torus inside within 
r = 35, while at t ~ 18,000, the mass accreted was equivalent to 
the torus mass out to r ^ 40. If accretion occurs radius after ra¬ 
dius, i.e., without any radial mixing, then gas located initially near 
these radii should fall on the BH during these periods of time. The 
initial Bernoulli numbers at r = 35 and r = 40 were Be = -0.012 
and -0.010, respectively. The meeasured total efficiencies averaged 
over the corresponding periods were y/totai = 0.029 and rj - 0.031. 
respectively. If the total efficiency of energy extraction was to re¬ 
flect the binding energy of the initial gas then it should drop by 
0.002 between the first and second periods (less bound gas accreted 
on the BH effecively deposits less energy at inflnity). However, the 
efficiency increased by a similar amount. We conclude that the gas 
mixes effectively before reaching the BH and that the measured to¬ 
tal luminosities are not sensitive to how much the initial torus was 
bound. 

The total luminosity indicates how much energy will be ulti¬ 
mately deposited at inflnity, where binding energy is zero and there¬ 
fore the outflowing flux of binding energy must have converted into 
other forms of energy. For supermassive black holes (SMBHs), we 
may expect that all this energy will ultimately affect the interstel¬ 
lar medium around the galactic nucleus and will contribute to AGN 
feedback. 

In analogy with equation (|^, we can straightforwardly de¬ 
fine an equivalent total radiative luminosity, i.e., radiative flux inte¬ 
grated over the whole sphere. 


n 2n 


( 11 ) 


which describes the radiative component of the total luminosity 
Ltotai- However, this quantity is not as fundamental as Ltotai since 
it is no longer conserved. In particular, radiation can be emitted 


or absorbed by the gas and can also gain/lose energy via momen¬ 
tum transfer to the gas. Nevertheless, for certain limited purposes, 
^rad,total cau be uscful. Radial profiles are shown with dashed lines 
in Fig.|^ For all the simulations, Lj-ad,total is negative in the inner re¬ 
gion, increases with radius, and ultimately becomes positive. Neg¬ 
ative values correspond to regions where more photons are dragged 
(advected) with the flow inward than manage to escape. These are 
the regions where the photon-trapping effect (discussed in detail in 
Section dominates. The fact that the total radiative luminosity 
becomes ultimately positive reflects the fact that the flow is slower 
at larger radii and so it is easier for photons to escape there. 


4.2 Optically thin and outflow regions 

Because of the limited range of inflow/outflow equilibrium in the 
disk mid-plane, extending at best only up to req ~ 25, it is im¬ 
possible to determine say how the outflowing energy is Anally dis¬ 
tributed between radiation and other forms when it reaches inflnity. 
However, the funnel and wind regions are converged to larger radii 
because of their higher velocities which allow them to be causaly 
connected with the equilibrium innermost region near the equato¬ 
rial plane. This fact allows us to measure luminosities in the funnel 
to larger distances from the BH. 

We divide the outflow region into two zones: (i) an optically 
thin zone which is visible to an observer at inflnity, and (ii) an out¬ 
flow region where the gas is on average energetic enough to reach 
inflnity. In all the simulations, zone (i) is a subset of zone (ii). The 
border of this zone is deflned to be the photosphere, which satisfies 
the following conditioij^ 



i.e., the total optical depth integrated along fixed polar angle from 
the domain boundary equals 2/3. The outflow region is deflned as 
the region where the relativistic Bernoulli parameter (eq.[^ is pos¬ 
itive. The borders of the two regions are denoted by dashed blue 
(optically thin) and green (outflow) lines, respectively, in Fig. 
Only for model rOOl (and its axisymmetric counterpart d300) does 
the optically thin region extend down to the BH. For all the other 
simulations the density of gas in the funnel region is large enough 
to move the lower edge of the photosphere away from the BH. For 
model r003, which is characterized by a significantly larger ac¬ 
cretion rate, the photosphere formally is at the outer edge of the 
domain. The photosphere is far from the BH also for model rOll 
for which the accretion rate (in Eddington units) is almost twice 
as high as in the fiducial one. The photosphere is relatively close 
(r ^ 40 at the axis) for model r020, which has a similar accretion 
rate as rOOl, but a higher BH mass. |McKinney et~n7| ( |2015| ) has re¬ 
cently showed that if strong magnetic held is present near the axis, 
one may get optically thin funnel down to the BH even for highly 
super-critical accretion rates. 

Radiation flowing out in the optically thin region is guaran¬ 
teed to reach observers at inflnity - no signiflcant interaction with 
gas is taking place here. The radiative luminosity integrated over 
this region may thus be viewed as a lower limit on the total radia¬ 
tive luminosity. The photons trapped in the optically thick regions 


^ This estimate of the optical depth is for a light ray propagating radially 
outward, where we have included the effect of the motion of the gas but 
neglected gravitational deflection of the ray. We also avoid edge effects, as 
described in |Sadowski et al.|j2015a| . 
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can in principle ultimately escape and could increase the radiative 
luminosity. However, this additional luminosity is not expected to 
exceed a couple of L^dd^ because if locally the radiative flux exceeds 
the Eddington flux L^ddlc^ and the gas is optically thick, this flux 
would accelerate the gas and convert radiative energy into kinetic 
energy of the outflo\\|^ 

To obtain the radiative luminosities of the optically thin and 
outflow zones, we calculate. 


-^rad, thin 
-^rad, out 



V^d0d6>, 

V^d0d6>, 


(13) 

(14) 


where ^thin and ^out denote the limits of the regions of optically thin 
and outflowing gas, respectively. The kinetic luminosities are cal¬ 
culated in a similar way 


Tkin, thin 
Tkin, out 


r thin 

I (u,+ yP^,)pu'yPgd(j)de, 
Jo 

r out 

I (Ut + 

Jo 


(15) 

(16) 


In the Newtonian limit one has: -(ut + ^|-gtt) l/2u^. 

Table |2] lists the radiative and kinetic luminosities in the two 
regions as measured at radii r = 50 and r - 250. For the flducial 
model rOOl (accreting at lO.OMsdd) only O.SLsdd escapes in the op¬ 
tically thin region at radius r = 50. This value increases to 0.95LEdd 
at radius r = 250 reflecting the fact that the optically thin region 
extends further from the axis and more radiation is able to enter the 
optically thin funnel. The radiative luminosity in the whole region 
of outflowing {Be > 0) gas is l.lSLEdd and 1.57LEdd at r = 50 and 
r = 250, respectively. At the same time kinetic luminosities are 
much smaller. Hardly any kinetic energy escapes in the optically 
thin region, mostly because of the negligible mass flux of outflow¬ 
ing gas there. In the whole outflow region, the kinetic luminosity 
grows from O.lVLEdd at radius r = 50 up to 0.35LEdd at r = 250. 
These numbers correspond to 13% and 18% of the total (radiative 
plus kinetic) luminosities, respectively. The increase of the frac¬ 
tional contribution of kinetic energy reflects the fact that radiation 
gradually accelerates gas, thereby losing momentum/energy. 

This effect is clear especially for simulation r003 (accreting at 
175MEdd)- The fractional contribution of kinetic luminosity grows 
from 39% to 58% between r = 50 and r = 250. The transfer of 
energy from radiation to gas is more effective because of higher 
optical depth in this run. For the same reason, there is no optically 
thin region in this simulation within the computational domain. The 
fact that both luminosities in simulation r003 grow significantly 
between the two radii arises from the fact that there is strong inflow 
of gas along the axis within r ^ 30 which prevents photons from 
escaping. Only further out on the axis does the standard, radiatively 
driven outflow region form. 

The kinetic and radiative luminosities measured in the thin and 
outflow regions for runs rOOl and r003 are significantly lower than 
the total luminosities measured according to equation 111 - At radius 
r = 250 and in the outflow region, these luminosities contribute 
only 37% and 10% of the total, respectively. Where does the rest of 
the luminosity go? It remains still in the optically thick gas in the 


^ In principle, since gas in the optically thick wind moves radially outward, 
beaming can allow a super-Eddington flux at the wind photosphere, as mea¬ 
sured in the lab frame. However, the velocities in our simulations are only 
a modest fraction of c, so the enhancement factor is unlikely to be large. 


disk interior, where the total energy budget is dominated by the out¬ 
flowing flux of magnetic (viscous stresses) and binding energy. The 
energy carried out in these channels (certainly, the binding energy) 
will be ultimately converted into radiative or kinetic energy by the 
time it reaches infinity. Unfortunately, the limited range of the equi¬ 
librium solution in our simulations prevents us from addressing this 
question directly. 

The total efficiency in simulation rOl 1, which is the only one 
with a rotating BH, is significantly higher than in the other sim¬ 
ulations. So are the radiative and kinetic luminosities integrated 
over the outflow region. The extra input of energy from the rotating 
BH makes the funnel region very energetic (see the bottom-middle 
panel of Fig. [^. At radius r - 250 the energy flux in the funnel and 
optically thick outflow is dominated by radiative luminosity equal 
to 8.4LEdd- Most of these photons are carried with the gas and ul¬ 
timately will escape when they reach the photosphere, which for 
most angles is outside the computational domain. However, the ra¬ 
diative luminosity may decrease if radiation keeps transferring its 
momentum to the gas. The kinetic component of the luminosity is 
significant already at this radius (r = 250) — roughly 1.5LEdd is 
carried in the form of outflow kinetic energy. 

The remaining three-dimensional run (r020), which differs 
from the flducial run (rOOl) in the mass of the BH (IOOOMq in¬ 
stead of lOM©) shows very similar properties in all respects. 


4.3 Angular distribution of energy fluxes 

Figure shows the angular distribution of the (time- and azimuth- 
averaged and symmetrized) radiative and kinetic fluxes of energy 
in the funnel/outflow region for three simulations (top to bottom): 
rOOl, rOll, and r003. Blue and red lines correspond to fluxes 
measured at r = 50 and r = 250, respectively. Solid and dashed 
lines denote the radiative and kinetic energies, respectively. 

For the flducial model (top panel, M - lOMEdd) the energy 
flux (as discussed in the previous section) is dominated by radia¬ 
tion. The angular distribution of radiative flux follows roughly a 
Gaussian with half-maximum width at 0 = 0.35rad. The maximal 
flux at radius r = 50 is F ^ 12FEdd; h increases to F ^ 19FEdd at 
r - 250. The radiative fluxes decline significantly with increasing 
polar angle 6. For an observer at 6 = 0.5rad (~ 30°) the observed 
flux (and the inferred source luminosity) is only ~ 4FEdd- 

The numbers given above are meaningful for the optically thin 
region but less so for other angles where the wind is optically thick. 
In these regions, there is still significant interaction between gas 
and radiation. The radiation is likely to accelerate the gas further 
and the radiative flux is expected to decrease towards the Eddington 
limit. To study this effect quantitatively one would have to perform 
simulations in a much bigger box and for a much longer time. 

At radius r = 50 there is hardly any kinetic luminosity in the 
polar region of the flducial run (rOOl). Only further from the BH 
is the gas is accelerated. In contrast to the radiative flux, the kinetic 
energy flux is not concentrated at the polar axis but rather in a shell 
around 6 - 0.35rad, sim ilar to the jet/wind boundary discussed in 
S^dowski et al. ( 2013b| ). For the flducial model accreting at lOMEdd 
the maximal kinetic flux at that inclination equals ~ 4FEdd- 

The angular profiles for the run with a rotating BH (rO 11) look 
qualitatively similar to the profiles of the flducial run. However, the 
magnitudes of the fluxes are much higher, reflecting the increased 
efficiency of accretion due to energy being extracted from the BH. 
The maximal radiative and kinetic fluxes at the axis at r = 250 
exceed lOOFEdd and 20FEdd. respectively. 

A larger mass accretion rate changes the picture significantly. 
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Figure 6. Top row of panels: Logarithm of average gas density (colors) and streamlines of average gas velocity (thickness indicates the magnitude of the 
velocity) for five models: rOOl (leftmost), r003, rOll, r020, and d300 (rightmost panel). Model r020 corresponds to Mbh = IGOOM©), and so its density 
was scaled up by a factor of 100 to enable direct comparison with the other simulations, which had Mbh = lOM©. The blue dashed contour reflects the location 
of the scattering photosphere as seen from infinity along fixed polar angle. The green dashed contour separates the bound {Be < 0) and unbound {Be > 0) gas. 
Bottom row of panels: Logarithm of the magnitude of radiation flux and its streamlines (thickness indicates the magnitude of the flux). For model r020, the 
flux was scaled up by 100. 


The third panel of Fig.[7]shows the angular profiles of fluxes for run 
r003. The radiative flux is still beamed at the axis, and the maximal 
flux exceeds 40FEdd at r = 250. The kinetic energy flux has a dif¬ 
ferent shape than it used to. It is no longer concentrated in a shell 
but increases with angle, reflecting the fact that the wind carries 
a significant amount of energy over a wide solid angle. However, 
even at the axis, the kinetic energy flux is much higher than in the 
fiducial model — it now exceeds lO^Edd at r = 250. The kinetic 
flux in the funnel is a result of radiative energy being converted 
into kinetic energy within the optically thick region dSadowski &| 

|Narayan|20l5bt - 

The angular distribution of energy fluxes discussed above 
should be considered only approximate due to the limitations of 


the Ml closure scheme adopted in this work. Most importantly, 
we evolve only the first moments of the radiation field instead of 
evolving specific intensities directly as in [Jiang et al.| ( |2014^ and 
jOhsuga & Takaha^ ( |2015j ). The radiation observed by a distant 
observer should, in principle, be calculated as an integral of the 
specific intensity pointing towards the observer over the whole ac¬ 
cretion disk. The local radial flux gives only an approximation of 
this quantity. Furthermore, the Ml closure is known to have dif¬ 
ficulties in treating the region closest to the polar axis. We have 
substantially mitigated this problem by including an extra radiative 
viscosity ( jS^dowski et al.|2015aj l, but the coefficients involved had 
to be chosen somewhat arbitrarily. 
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Table 2. Luminosities 




rOOl 



r003 


rOll 



r020 


d300 (2D) 

Ml MEdd 


10.0 



175.8 


17.6 



10.1 


8.9 

f^total / I^Edd 


5.26 



83.27 


14.27 



5.14 


5.31 

7total 


0.030 



0.027 


0.084 



0.029 


0.034 

Hum 

50 


250 

50 

250 

50 


250 

50 


250 

50 250 

T^rad, thin / f^Edd 

0.30 


0.95 

0.0 

0.0 

0.0 


2.13 

0.20 


0.85 

0.64 2.14 

I^kin, thin I^Edd 

0.00 


0.15 

0.0 

0.0 

0.0 


0.16 

0.00 


0.15 

0.01 0.23 

f^rad, out. / I^Edd 

1.13 


1.57 

0.78 

3.47 

6.70 


8.39 

1.11 


1.60 

1.86 2.58 

I^kin,out/I^Edd 

0.17 


0.35 

0.50 

4.82 

1.41 


1.54 

0.21 


0.39 

0.31 0.41 


M - average accretion rate 

Ltotai/^Edd - total (conserved) luminosity of the system including the flux of binding energy 
^totai = Aotai/^c^ - efficiency of accretion in total luminosity 

^rad,thin? f-kin,thin ” radiative and kinetic luminosities in the optically thin polar region integrated at radius rium 
^rad,out 5 ^kin,out" radiative and kinetic luminosities in the region of unbound gas (Be > 0) integrated at radius rium 



^3 






1.0 


Figure 7. Average energy fluxes in simulations rOOl (top panel), r003 
(middle), and rOll (bottom) as a function of the polar angle 6. Blue and 
red fines correspond to fluxes measured at r = 50 and r = 250, respectively. 
Solid fines show the radiative flux and dashed fines show the flux of kinetic 
energy. Bright fine segments are within the outflow region (Be > 0), while 
the dull shaded segments are outside this region. 


5 PHOTON TRAPPING 

At the time the thin disk theory was established ( [Shakura & Sun-| 
|yaev|197^|Novikov & Thorne|1973| l, it was commonly understood 
that the standard radiation pressure dominated disk can extend up to 
and above the critical Eddington accretion rate (equation [^. Once 
the Eddington limit is exceeded, it was predicted that the most ener¬ 
getic, inner region would attempt to produce a luminosity exceed¬ 
ing the Eddington value and this would produce a radiatively driven 
outflow inside the so called spherization radius. At the same time 
the wind would modify the accretion rate on the BH. 

This picture did not take an important effect into account. 
When the accretion flow is very optically thick, photons do not have 
enough time to diffuse vertically to the disk photosphere before 
they are dragged radially inward by the accreting gas and advected 
into the BH. This effect was described in a simple spherical context 


by|Begelman|(|1978 

1 and included for the first time in a full-fledged 

accretion model by 

Abramowicz et aL|(|1988|) who constructed the 

so-called slim disk model. These authors predicted that the radia- 


tive luminosity would no longer scale with the accretion rate above 
the critical rate. Also, in priniciple, the slim disk state would pre¬ 
vent the spherization phenomenon. 

As we have shown in the previous section, and as we will ex¬ 
plain in greater detail here, the truth is in between — there is a 
region, but only near the axis, where locally flux is significantly 
super-Eddington and where radiation drives gas out of the disk, but 
at the same time photons are trapped and advected towards the BH 
in the inflow region near the equatorial plane. 

In this section we try to answer the question: where is photon 
trapping effective or, in other words, where is the trapping radius, 
the border between the radiatively inefficient (slim disk) and radia¬ 
tively efficient (thin disk) regimes? 


5.1 One-dimensional, luminosity-based trapping radius 

It is relatively straightforward to define the trapping radius in a 
spherically symmetric flow ( |Begelman|I978| ). Assuming a Bondi- 
like flow, we can compare the local outward radiative diffusion ve¬ 
locity to the local gas inflow velocity and thereby estmate 

-^trap, Bondi ~ ^I ^Edd- (l'^) 

This approach gives only an upper limit on the location of the trap¬ 
ping radius in a real accretion flow; accretion disks are not spheri¬ 
cally symmetric but have low density polar regions where radiation 
can more easily escape, and also the inflow velocity near the equa- 
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torial plane is significantly lower than the free-fall velocity assumed 
in the Bondi model ( [Sqdowski et al.|2014] ). 

The simplest way of estimating the location of the trapping 
radius numerically is to look at the net fiux of radiation, e.g., the 
dashed lines in Fig. Negative values mean that more photons 
were moving inward (mostly because they are dragged by optically 
thick gas) than outward. The radius at which the radiative lumi¬ 
nosity goes to zero could then be defined as the trapping radius. 
This approach is simple and provides a single effective trapping 
radius. It does not, however, account for the non-uniform struc¬ 
ture of the disk - photons easily leave the system at the axis and 
are more effectively trapped near the disk mid-plane. The present 
simple definition weights both effects and provides some kind of 
an average. By this definition, for simulations rOOl and r020, the 
effective trapping radius is located around r - 35. Model r003, 
which has a significantly larger accretion rate, has almost the same 
trapping radius, r 40. Model rOll, with a rotating BH, produces 
much more powerful radiation fiux along the axis, and correspond¬ 
ingly the effective trapping radius is much closer to the BH (r ^ 10, 
although the trapping in the bulk of the disk is as effective as in the 
other cases). A caveat: The values given above for the simulations 
with non-rotating BHs should be taken with caution, because the 
bulk of the disk is outside the infiow/equlibrium region at these 
radii. 


5.2 Two-dimensional trapping - importance of the diffusive 
flux 


Because of the limitations of the definition given in the previous 
section, we discuss here a different approach which, in particular, 
seeks to account for local properties of the flow. 

I Jiang et ar] ( |2014b| ) calculated energy-weighted, vertically in¬ 
tegrated, radial and vertical velocities of radiation transport. By 
comparing the two one can distinguish whether more energy flows 
up (away from the midplane) or inward (towards the BH). not ade¬ 
quate to say where the photon trapping is effective because it does 
not gives no importance to diffusion. The vertical flow of radiation 
could be the result of advection (in a wind) or turbulence | Jiang et| 
[n^ ( |2014b] ). Thus, while it is a good way of comparing vertical and 
radial energy transfer, it provides no information on the transport 
mechanism. 

Below we attempt to quantify photon trapping by calculating 
the fraction of the total radiative fiux that comes from diffusive 
transfer. To estimate the diffusive flux we use the moment equa¬ 
tions 1^, assuming dt - 0 and a dia^nal form of the radiative 
stress energy tensor in the fluid frame (^" = l/3^^0» to obtain 




1 


3(/^'p) dx^ 


AE), 


(18) 


which is the standard diffusive flux formula ( [Rybicki & Lightm^ 
|1979| ). This expression should be used only inside the optically 
thick parts of the disk where the Rosseland approximation is sat¬ 
isfied. 

The top panel of Fig. shows the magnitude (colors) and di¬ 
rection (red streamlines) of the total radiative fiux, = RJ, for the 
fiducial model rOOl which accretes at lOMEdd- Deep in the disk the 
radiation flows towards the BH, while the polar region is filled with 
optically thin radiation escaping along the axis. There is a transi¬ 
tion region at intermediate angles, 6 ^ 45°, where radiation rela¬ 
tively smoothly switches from flowing radially inward to flowing 
outward. 
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Figure 8. Top panel: Magnitude and streamlines of the total radiative flux 
F'ad = -FJ. Middle pan el: Magnitude and direction of the radiative diffusive 
flux (equation II^. Bottom panel: Magnitude of the comoving frame 
radiative energy density (colors) and the ratio of the diffusive radiative flux 
to the total radiative flux (F^-g:/F*a^). All panels correspond to simulation 
rOOl (lOMEdd)- The yellow stars denote the location where we studied in 
detail the vertical fluxes plotted in Fig.[7o| 
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Figure 9. Magnitude of the comoving frame radiative energy density (col¬ 
ors) and the ratio of the diffusive to advective radiation flux 
simulation r003 (ITSMEdd)- 


Blue streamlines in the same panel show the average gas ve¬ 
locity. If all the radiative flux was coming from photons advected 
with the gas, the total radiation flux vectors would follow every¬ 
where the direction of the gas velocity. The agreement between the 
two sets of streamlines is very deep inside the disk, where both gas 
and radiation flow inward. This is where gas is most optically thick 
and where one expects efficient photon trapping. The streamlines 
agree again in the polar region, but this is not because of photon 
trapping. The low optical depth allows radiation to flow indepen¬ 
dently of the gas. The funnel geometry, however, makes both gas 
and radiation to flow upward, and the locally super-Eddington flux 
accelerates gas in its direction. In the intermediate region between 
the funnel and the disk interior, the gas velocity and radiation flux 
vectors do not point in the same direction. As mentioned previ¬ 
ously, gas flows on average radially inward in the disk and only 
close to the surface is it blown away, causing the velocity stream¬ 
lines to turn rapidly outward. Radiation flux, on the other hand, 
changes direction rather smoothly. This difference between gas and 
radiation streamlins shows that there must be an additional compo¬ 
nent of the total flux besides the advective photon transport. This is 
of course diffusive transport. 

The middle panel in Fig. shows on the same axes and with 
the same color scale the diffusive flux estimated according to equa¬ 
tion GD- As expected, the diffusive flux follows the gradient of 
radiative energy density, i.e., radiation diffuses in the vertical di¬ 
rection away from the equatorial plane. Then, it turns smoothly to¬ 
wards the axis and enters the funnel region. IN the funnel itself (the 
red region), the estimate of the diffusive flux must be disregarded 
since gas there is optically thin and radiation is free streaming in¬ 
stead of diffusing. However, the diffusive flux estimates are trust¬ 
worthy in the disk interior and the transition region between the 
disk and the funnel. 

We now check whether the diffusive flux is, in fact, responsi¬ 
ble for the deviation between the total radiative flux streamlines and 
the advected radiation flux streamlines (which would be aligned 
with the gas velocity). For this purpose we choose a location in the 
region of largest deviation, {r,0) - (20,51°) (denoted by the yel¬ 
low stars in Fig. |^, and study local fluxes of radiation over time. 
Fig. [^presents the estimated orthonormal polar component of the 



Figure 10. Estimated polar diffusive flux of radiation (equation [T^ ver¬ 
sus the excess of the total flux over the local advective flux of radiation 
(equation]^. Colors denote the local optical depth over a distance of one 
gravitational radius. The dashed line shows where the two fluxes agree, i.e., 
when the excess can be explained entirely by diffusive transport. The data 
are from simulation rOOl and the points range from t - 1 ,500 20,000 

with a cadence of dt = 50. 


diffusive flux, calculated as. 


1 


— _ j7 

^diff “ O_ 


(19) 


3kp rdO 

as a function of the difference between the polar components of the 
total flux and the estimated advective flux. 


po _ pO - /?< 


-Eu^\r. 


( 20 ) 


The signs have been chosen such that the positive values corre¬ 
spond to fluxes pointed out of the disk, i.e., towards the axis. The 
sizes of markers measure the local optical depth for scattering as es¬ 
timated by Tes = KqsPRg- The dashed line denotes = T’^ad“^fdv- 
The plotted points cover the time period t = 7,500 ^ 20,000 with 
cadence of At = 50. 

The majority of points in Fig. cluster around the dashed 
line, showing that the excess of radiative flux over advective flux 
has exactly the analytically predicted magnitude, i.e., the excess of 
polar radiative flux is due to diffusion. The agreement is no longer 
perfect at the lowest optical depths. This is natural since, when the 
optical depth is low, radiation decouples from the gas and the diffu¬ 
sive approximation no longer valid. In fact, this is where we expect 
the actual flux to lie below the analytically estimated diffusive flux, 
as is indeed seen in the plot. 

The colors of the markers in Fig. [^denote the ratio of mag¬ 
nitudes of instantenous diffusive and advective fluxes. Despite the 
fact that the time-averaged advective flux has hardly any polar com¬ 
ponent, instantenous advective flux in 6 direction can be more than 
30 times stronger than the corresponding diffusive flux. However, 
because of turbulent motion of gas, it averages to a value signifi¬ 
cantly smaller than the diffusive flux. Therefore, it is the diffusive 
flux which dominates the average net radiation flux towards the axis 
near the surface of the disk. 

The bottom panel in Fig. shows the ratio of the magnitude 
of the diffusive flux over the total flux plotted with contours on top 
of the distribution of fluid frame radiative energy density, E. Deep 
in the disk, advection dominates strongly over diffusion. The mag¬ 
nitude of the former flux is > 30 times larger than the magnitude 
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of the diffusive flux. This ratio is lower closer to the disk surface, 
where the density drops and the gas is no longer able to drag pho¬ 
tons efficiently. In the transition region discussed above, the ad- 
vective and diffusive fluxes have comparable magnitudes, with the 
diffusive flux dominating the polar component. 

Figure shows the same ratio of the diffusive to advective 
fluxes, but for simulation r003 which has almost 20 times larger 
accretion rate. The larger accretion rate implies higher gas density 
and optical depth, and, as a result, more effective photon trapping. 
This is indeed the case. The advective radiative flux is now > 300 
times larger than diffusive in the disk near the equatorial plane. 

The properties described above show that for moderately 
super-critical accretion rates, ~ lOMEdd. photon transport deep in¬ 
side the disk is dominated by radial advection. However, diffusive 
transport of energy becomes important further from the equatorial 
plane, near the transition between the disk and the funnel, where the 
gradient of the energy density is large, and the optical depth of the 
gas is no longer huge. Radiation in the optically thin funnel follows 
the axis and is decoupled from gas. Because of significant photon 
trapping deep in the disk, only photons generated close to the sur¬ 
face can escape the disk and join the funnel region. This explains 
the relatively low radiative luminosities of our simulated disks and 
is in agreement with the slim disk model (e.g., |Abramowicz et al.| 
|1988[ |S^dowski||2011| |. For higher accretion rates, the region of 
dominant photon trapping extends further out, not only near the 
equatorial plane, but also in the (now optically thick) funnel. Like¬ 
wise, one expects that, for lower accretion rates, photon trapping 
will become less and less effective and the radiative efficiency will 
approach the thin disk value. 


5.3 Turbulent transport 


So far we have shown that the diffusive transport of radiation is 
important near the surface of the disk, and that the radial photon 
advection dominates over diffusion deep inside the disk. There is 
potentially one more way of radiation transfer - turbulent transport 
pointed out by |Jiang et al.| ( |2014b] ). It is effective if radiative energy 
is transported without transporting mass, similar to what happens 
in convection. Such behavior may result, e.g., from magnetic buoy¬ 
ancy. To assess the importance of this effect it is enough to compare 
the density- and radiative energy-weighted vertical velocitities of 
the gas. If the former is zero and the latter is not, then radiation is 
transported by gas without transporting mass, and hence, it is dif¬ 
ferent in nature from the standard photon trapping which results 
from mean motion of gas. 

If turbulent transport is effective then lighter gas, but contain¬ 
ing more radiation, moves preferably towards the disk surface. In 
Fig. we plot the polar velocity of gas as a function of radiation 
to rest-mass energy density ratio for the same location as in Fig.p^ 
(red), and for another point closer to the equatorial plane, located 
at (r, 6) = (20,80°) (blue markers). The size of the markers denotes 
the optical depth as it did in Fig.[^ We see that there is no corre¬ 
lation between the relative radiative energy density content and the 
vertical motion of the gas, what suggests that the turbulent effect 
cannot be strong. 

From the same sets of points we now calculate the density- 
and energy-weighted velocities, defined as. 
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Figure 11. Vertical velocity of gas as a function of radiative to rest-mass en¬ 
ergy ratio for simulation rOOl and two points, located at (r, 9) = (20, 851°) 
(red) and (r, 6) = (20,80°) (blue markers). 


where the sums go through all the points in the set, corresponding 
to different moments of time and fixed location. For the point closer 
to the disk surface we get {v^)p - 0.0022 and = 0.0013. For 
the other point, located almost at the equatorial plane, we have, 
{v^)p = 0.0038 and = 0.0002. All the values are positive 
what indicates that both the gas and radiation is advected with the 
gas towards the equatorial plane. However, the magnitudes of these 
velocities are at least order of magnitude lower than of the cor¬ 
responding radial velocities, what reflects the fact that the radial 
motion and advection dominates. One may, however, notice that 
the energy-weighted velocities are lower, i.e., they tend to deviate 
from the density-weighted velocities towards the surface. Never¬ 
theless, the magnitude of the turbulent transport of radiation is not 
significant. 


5.4 Advection coefficient 

Defining the efficiency of photon trapping in a non-spherical ac¬ 
cretion flow is a challenging task because of the multiple dimen¬ 
sions involved. Here we try to calculate the advection coefficient 
which estimates the fraction of photons generated at given radius 
that are able to escape, the rest being advected to the BH. For this 
purpuse, at each radius rbox, we consider a disk annulus extending 
from ri = rbox - O.Svq to r 2 = rbox + O.Sro, and limited in 6 to the 
range 6+ - 90° + 37° (Fig.[^. The polar angle range is chosen to fit 
the center of the transition region discussed in the previous section 
where the radiation flux in the fiducial model is dominated by its 
polar component. The annulus covers extends over In in azimuthal 
angle. From the time averaged disk properties we extract luminosi¬ 
ties of the radiative flux that crosses each surface of the annulus. 
For the radial sides, we compute 

^6+ 

R]^rkd0d4,, (23) 

0 Jo- 

where the integration takes place either at the inner edge r - ri ox 
at the outer edge r - r 2 . The fluxes crossing the top and bottom 
surfaces are similarly integrated to give, 

J f*2n ^r2 

I R1^Pgdrd(p, (24) 

0 Jri 
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Figure 12. Shows the box and the definitions of the luminosities used in 

equation 0 to estimate the advective factor ^adv = -f-r, 2 + 1 - 0 ) • 


where the integration is done at fixed polar angles, 6+. 

The integrated fiux L^ i incoming through the outer radial 
sufrace tells how much radiation is advected into the volume. The 
corresponding luminosity Lr^i crossing the inner surface tells how 
much radiation is advected out of the same volume. Meanwhile, the 
total radiation generated within the annulus is i + Lg- L^ 2 - Thus 
Lr^i - L ^2 divided by the latter quantity measures the fraction of 
energy that is advected radially, whereas Lq over the same quantity 
measures what fraction of the radiative energy escapes through the 
top and bottom surfaces. This then motivates the following defini¬ 
tion of the advection coefficient, ^adv. 


^adv ~ 


Lr,\ - Lr^2 
Lr,\ - Lr^l + Lq 


(25) 


In the limit of a radiatively efficient disk we would have = 
L ^2 ~ 0 and ^adv ~ 0- In the opposite limit of an advection- 
dominated flow we expect Lq = 0 and ^adv ~ 1- Tehrefore, it is 
natural to define the effective trapping radius as the location where 
^adv = 1/2, i.e., where half the radiation generated at that radius 
manages to escape and half is advected to the BH. 

Figurep^shows the above advection coefficient ^adv as a func¬ 
tion of radius for the fiducial model rOO 1 and the large accretion 
rate model r003 (the chosen vertical size of the box corresponds to 
the region of purely polar radiative flux only for the former). For 
the fiducial simulation, the advection factor increases towards the 
BH, as expected because of increasing inflow velocity. At radius 
r = 10 we find ^adv = 0.8 which means that only ~ 20% of pho¬ 
tons generated at that radius manage to escape and enter the funnel. 
The coefficient drops down to ^adv = 0-55 at the edge of the in¬ 
flow/outflow equilibrium region (r = 25). The profile suggests that 
the effective photon trapping radius, defined by ^adv = 1 /2, is prob¬ 
ably around r ^ 35 (the point plotted at r = 30 is outside the region 
of inflow equilibrium and is a little suspect). 

The efficiency of advection for the simulation with higher ac¬ 
cretion rate (r003) is significantly larger because of the larger op¬ 
tical depth. Even at radius r - 25 only ~ 15% of photons escape 
the disk. In the innermost regions of this model, ^adv exceeds 1, re¬ 
flecting the fact that the top-bottom luminosities are negative, i.e., 
photons are brought into the box, and no radiation escapes. 

In our analysis in this subsection we focused solely on the 



R/M 

Figure 13. The advection coefficient, ^adv (equation as a function of 
radius for simulations rOOl (blue) and r003 (red circles). 


radiative energy flux. In general, one should include also the flux of 
mechanical energy because dissipation may result in kinetic energy 
leaving the box. Including this component hardly affects the values 
of the advection coefficient calculated above. 


6 VARIABILITY 

Radiation coming from accretion disks is known to be highly vari¬ 
able (e.g., |Done et al.|2007| l. In case of galactic BH binaries, this 
variability takes place on short timescales (the horizon crossing 
time for a lOM© BH is GMjc^ ^ 5 X 10“^s). The variability is 
strongest in the hard state, and weakest in the thermal state. How¬ 
ever, even in the latter, the power spectrum is far from featureless. 

Studying variability is a powerful tool in understanding accre¬ 
tion flows. The characteristic frequencies tell us where the modu¬ 
lated radiation come from. Features in the power-spectrum, e.g., 
quasi-periodic oscillations or breaks, can manifest more subtle 
properties of the disks (e.g. |Ingram et al.|2()0^|Wellons et al.|20I4| |. 
Modeling the variability and its power spectrum has so far been 
limited mostly to analytical models. However, the dynamics of the 
gas and the properties of the radiation field are complicated and 
highly non-stationary, so the analytical approach is limited, and ul¬ 
timately we should model variability using time-dependent three- 
dimensional simulations. 

However, this approach is not straightforward. Numerical 
modeling of radiation in MHD codes is limited by a number of 
factors. First of all, radiation is solved for in the grey approxima¬ 
tion and some arbitrary (usually blackbody) shape is often assumed 
for the spectrum. Secondly, general relativistic effects are rarely in¬ 
cluded. Moreover, Comptonization is either neglected or treated in 
a crude way. Last but not least, various approximations for radiation 
closure are adopted (with the exception of direct radiation transfer 
solvers operating on a fixed grid of angles, as implemented recently 
by |Jiang et al.| ( |2014^ and |Ohsuga & Takahashi| ( |20I5] )). Simplistic 
closure schemes may limit the information available for calculating 
the visible spectra. 

The most reasonable way to proceed is to take the global, time- 
dependent output of a disk simulation and postprocess it with a so¬ 
phisticated radiation (and only radiation) solver which will not be 
as limited as full radiation-MHD simulations. Such codes, which 
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Table 3. Fractional variability of radiative flux in model rOOl 



6 = 2° 

6 = 21° 

o 

o 

II 

0.20 

1.96 

r = 250 

0.23 

0.64 

r = 500 

0.23 

0.58 


solving the frequency-dependent radiative transfer equation and ac¬ 
count for relativistic effects and Comptonization have recently been 
developed (e.g. |Zhu et al.|2015] l and are expected to be soon avail¬ 
able for spectral modeling. 

In the meantime, we attempt to directly estimate the variabil¬ 
ity of light curves from the simulations with KORAL described in this 
paper. Because of the limitations mentioned above, in particular the 
fact that we evolve only the first moments of the radiation (Ml clo¬ 
sure), our approach is expected to provide only a rough qualitative 
understanding of the temporal properties of radiation coming from 
super-Eddington accretion flows. 


6.1 Light curves 

In principle, the observer is located at infinity, i.e, r » rout- The 
correct way of measuring the radiation reaching a distant observer 
is to integrate the specific intensities pointing towards the observer 
at the outer edge of the computational box. However, because of 
the limited range of inflow/outfiow equilibrium near the equatorial 
plane, we are limited only to studying light escaping in the polar 
region near the axis, where the disk solution has converged to large 
enough radii. Moreover, the size of our computational box is ob¬ 
viously limited and we cannot measure the radiative fluxes at radii 
larger than r > rout/2 = 500. Even if the duration of simulation 
was infinite, the photosphere would be located beyond the domain 
boundary except in the polar region, and even there we resolve the 
photosphere only at the lowest accretion rates. 

Because of the limitations mentioned above, we decided to 
estimate the light curve by looking at the local flux of radiation 
at some location inside the polar region. To be sure that the light 
curve measured in this way is close to the light curve seen from in¬ 
finity, we insist that the radiation is already decoupled from the gas; 
otherwise, continued interaction with (and acceleration of) the gas 
could drastically change the properties of the radiation that finally 
escapes. Once this condition is satisfied, and provided there is little 
radiation coming from other regions of the disk towards the chosen 
line of sight, the shape of the light curve should be independent of 
the radius it is measured at, only shifted in time by the light cross¬ 
ing time. Below we test if this criterion is satisfied for the locally 
measured fluxes in simulation rOOl. 

To extract the light curves we choose an arbitrary slice through 
the poloidal plane, define the inclination angle, and choose a radius 
where we measure the energy fluxes. Then we take the radial com¬ 
ponent of the lab-frame radiative flux, and plot it as a function 
of time. Eig.p^ shows such time series extracted at two polar an¬ 
gles {6 - 2° and 27°, top and bottom panels, respectively), and at 
three radii (r = 100, 250, and 500). Only the second half of the 
simulation, when the accretion has settled down to a quasi-steady 
state (compare Eig.[^, is shown. The fluxes measured at r = 100 
and r - 500 were scaled to account for geometrical expansion to 
match the magnitude of the flux at r = 250. 

The first panel corresponds to radiation escaping along the 
axis. Eor model rOOl, gas is optically thin on the axis at all radii 
so one may expect the measured fluxes not to depend on the ra¬ 
dius where the measument is performed. This is indeed the case. 



Figure 14. Variability of the radiative flux in model rOO 1 measured at polar 
angle 6 = 2° (top panel) and 21° (bottom panel) and at three radii: r = 100, 
250, and 500. 


The profiles resemble each other to a good accuracy, at least for 
t > 12000, with the shift in time reflecting the propagation of light 
in the vertical direction. This gives some confidence that the light 
curve we calculate is what a distant observer might see. It is in¬ 
teresting, however, that the propagation speed down the funnel to 
match the time delay between radii is less then the speed of light. It 
equals approximately 0.25c which is the characteristic velocity of 
the radiation field (the velocity of the radiation rest frame) inside 
the funnel. It is determined by the geometry of the funnel and the 
effect was discussed close to 40 years ago by |Sikora] ( |1981| ). Ba¬ 
sically, not all the photons go straight up; a significant fraction go 
sideways and some even go backward. The characteristic velocity 
is in some sense the average radial photon velocity. 

The short duration of the light curves (~ 5^ real time for a 
IOMq BH)) do not allow us to calculate power spectra. Instead, we 
calculate the fractional variability through 


mm ’ 


(26) 


where cr denotes the standard deviation and () stands for the mean 
luminosity. The second column of Tablej^gives the fractional vari¬ 
ability calculated for light curves in the top panel of Eig. Note 
that the fractional variability does not change much with radius. In 
particular, it does not change between r = 250 and r = 500 which 
again reflects the fact that the radiation flowing along the axis has 
established its temporal profile and is not affected any more by in¬ 
teraction with gas or radiation at larger radii. 

The bottom panel of the same figure presents radiative light 
curves measured at the same radii but at a larger inclination an¬ 
gle, 6 = 27°, which corresponds to the edge of the funnel region, 
but located in the optically thick portion of the gas. One may ex¬ 
pect therefore, that the radiation interacts (is absorbed or emitted) 
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with the gas up to large radii, even out to the domain boundary, 
and that the radiation fluxes measured at different radii will differ, 
even after taking the geometrical factor into account. The radiation 
flux measured at r = 100 (blue line) varies significantly, reach¬ 
ing at moments negative values. This behavior reflects the motion 
of the gas to which the radiation is coupled to. The gas on aver¬ 
age moves out, but temporarily can move inward, dragging photons 
with it. The light curves become more smooth and the luminosities 
become larger in magnitude with increasing radius. These result 
from the fact that the gas moves in a more laminar way further out, 
and that radiation originating from a larger volume contributes to 
the measured luminosity. As a result, the corresponding fractional 
variability (third column of Table decreases significantly with 
increasing radius. 

Because of the significant coupling with gas, and the fact that 
we do not resolve the photosphere (because the computation box 
is not large enough box and the duration of the simulation is too 
short), the light curves for the larger inclination angle are not reli¬ 
able. Only the light curve for radiation leaving the system almost 
exactly along the axis is robust. However, even this will not be the 
case if the accretion rate is larger since then, even on the axis, the 
photosphere will be located at large radii. This is a severe limitation 
for studies of variability in super-critical disks. In contrast, for thin 
disks, the photosphere is close to the equatorial plane and one can 
hope to extract useful variability information from simulations. 


7 AXISYMMETRIC SIMULATIONS 

Angular momentum in accretion disks is transported by MHD tur¬ 
bulence, which is driven primarily by the magnetorotational insta¬ 
bility palbus & Hawley|1991| l. Turbulent motion in the poloidal 
plane leads not only to the exchange of angular momentum, but 
also to dissipation of magnetic held. In axisymmetric simulations 
this process quickly causes the magnetic held to decay, and at some 
point, shortly after the beginning of a simulation, typically after t ~ 
5000, the accretion stops. In reality, the poloidal magnetic held is 
revived by a dynamo process which results from breaking axisym- 
metry — non-axisymmetric perturbation can “rotate” toroidal held 
into the poloidal plane and thereby regenerate poloidal held. Three 
dimensional simulations of turbulent magneto-fluids in shearing 
boxes (and also in global models) have shown that the saturated 
state is characterized by an average toroidal-to-radial magnetic held 
ratio 6b ~ 0.25, and an average magnetic-to-gas pressure ratio 
p = 0.1. There is no way of obtaining such a saturated state in 
pure axisymmetric simulations. 

On the other hand, the advantages of assuming axisymmetry 
and simulating accretion disks in two dimensions are obvious - 
such simulations are cheaper by more than an order of magnitude in 
terms of the required computational time. Until recently, however, 
axisymmetric simulations were limited to extremely short durations 
because of the rapid decay of the magnetic held mentioned above. 
The situation has now changed. [S^dowski et al.| ( [2015a| ) introduced 
a mean-field model of the dynamo effect which mimics the prop¬ 
erties of three-dimensional MRI-driven turbulence but can be ap¬ 
plied to axisymmetric simulations. In this approach, the properties 
of the magnetic held are driven towards the prescribed characteris¬ 
tics of the saturated state, described by parameters 6b and p'. This 
is achieved by “pumping” new vector potential into the MHD flow, 
leading to a correction to the existing poloidal held. The poloidal 
held is enhanced in regions where the magnetic held is too weak, 
and the toroidal component of the held is damped in regions were 


the magnetic pressure exceeds the prescribed saturation value. It 
has been shown that such an approach successfully allows for arbi¬ 
trarily long simulations and indeed leads to a saturated state similar 
to that seen in three-dimensional simulations. 

Our aim in this section is to compare the properties of three- 
dimensional simulated super-critical, radiative disks with corre¬ 
sponding two-dimensional axisymmetric simulations. For this pur- 
puse we simulated an additional model (d300) which was run in 
axisymmetry, with 252x234 grid points in the poloidal plane, and 
initialized in an identical manner to the fiducial model rOO 1 (see 
Table [T). We used the mean-field dynamo model with the fiducial 
saturated state parameters 6b - 0.25 andyS' = 0.1. Model d300 was 
run until t - 190,000, an order of magnitude longer than the three- 
dimensional model rOO 1, yet it required less computer resources by 
almost an order of magnitude. This enormous saving is the reason 
why we feel it is important to explore the two-dimensional option 
fully. 

The left- and right-most panels in Fig. compare the den¬ 
sity distribution on the poloidal plane, gas velocities, radiative flux, 
location of the photosphere, and border of the outflow region in 
the three- (rOOl) and two-dimensional (d300) simulations. Quali¬ 
tatively, the two solutions (after averaging over time and azimuth 
in the three-dimensional case) agree very well. The gas shows the 
same dynamical properties, with outflow taking place only in the 
funnel region, the magnitude of the radiation flux is similar, and 
the photosphere is located at the same place. The only noticeable 
difference is near the disk surface at larger radii (r > 20) where 
the two-dimensional run shows more significant vertical motion. 
This may be because of the approximate treatment of the dynamo 
effect which is constructed to satisfy only the vertically averaged 
criteria, without paying too much attention to the vertical profile of 
magnetic held properties. It should also be kept in mind that the 
three-dimensional model has achieved inflow/outflow equilibrium 
only for r < 20. It could be that even the relatively minor differ¬ 
ences in structure between the three- and two-dimensional models 
would be reduced once the former is run for a much longer time and 
reaches equilibrium out to radii ~ 100. Very expensive simulations 
would be required to verify if this explanation is correct. 

To quantify the comparison between models rOOl and d300, 
we have calculated radial profiles of the accretion rate, surface den¬ 
sity, and density-weighted angular momentum for the two runs. 
They are shown respectively in the top, middle, and bottom panels 
in Fig.[^ The net accretion rates are similar. The two-dimensional 
simulation has M - 8.9MEdd^ roughly 10% lower than the three- 
dimensional model, but close enough to allow a meaningful com¬ 
parison of the two models. The extent of the flat section of the 
net accretion rate profile is a benchmark for the extent of in¬ 
flow/outflow equilibrium. As discussed earlier, for simulation rOO 1 
the equilibrium region extends only to r ~ 20. Thanks to the long 
duration of the d300 simulation, the equilibrium region here ex¬ 
tends much further out, to r ~ 100. Dashed lines in the same panel 
show the rate of mass lost in the wind. It is comparable for the two 
runs, and suggests that the amount of outflowing gas equals the net 
accretion rate around radius r - 60. 

The middle panel in Fig. [^compares the surface density pro¬ 
files. They are indistinguishable where the equilibrium regions of 
the two simulations overlap (the sections of the curves outside the 
equilibria are marked with shaded colors). Similarly, the profiles of 
angular momentum also match well (bottom panel) except in the 
plunging region inside the marginally stable orbit, where the two- 
dimensional run results in a slightly lower angular momentum. This 
difference probably arises from the fact that the dynamics of the 
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Figure 15. Radial profiles of the net and outfiowing accretion rate (top 
panel), surface density (middle), and specific angular momentum (bottom) 
for the three-dimensional run rOO 1 and the equivalent two-dimensional run 
d300. The light shaded line segments in the second and third panels corre¬ 
spond to regions where infiow/outfiow equilibrum has not been reached. 


flow and the properties of the magnetic held in this region are not 
governed by the combined effects of shear and the dynamo process 
( |Penna et al.|2013b| l. 

In summary, the above comparisons indicate that the average 
properties of the two- and three-dimensional simulations are re¬ 
markably similar. The temporal properties, on the other hand, are 
significantly different. In Fig. we show radiative light curves 
measured at angle 6-2° and at radius r = 250 (the light curve 
of simulation rOOl corresponds to the one shown in the top panel 
of Fig. for this particular radius). The two-dimensional run re¬ 



time [GM/c^] 


Figure 16. Time variability of the radiative flux measured at angle 6 = 2° 
from the axis at radius r = 250 for the three-dimensional simulation rOOl 
and the two-dimensional model d300. 

suited in noticeably larger radiative luminosity in the outflow re¬ 
gion and larger beaming towards the axis (but the same total lu¬ 
minosity, see Table [^. As a result, the flux in the axisymmetrical 
model is twice as large, and at the same time, it is also more vari¬ 
able - its fractional variability (/ = 0.64) significantly exceeds that 
of the three-dimensional run (/ = 0.23). The difference must come 
from the fact that the axisymmetrical simulation neglects the many 
non-symmetrical modes which develop in three-dimensional simu¬ 
lations. Variations in independent patches in azimuth would tend to 
wash each other out when averaged, but no such averaging would 
be present in an axisymmetric model. 


8 SUMMARY 

In this paper we presented a set of four three-dimensional simula¬ 
tions of black hole super-critical accretion disks. Two of the sim¬ 
ulations accreted roughly at lOMEdd (models rOOl and r020), the 
third simulation (r003) had a significantly larger accretion rate of 
176MEdd. and the fourth simulation (rOll) was characterized by a 
rotating BH with spin a* = 0.7 and accreted roughly at 17MEdd- 
The fifth simulation we presented was performed assuming ax- 
isymmetry (hence this was a two-dimensional simulation) and cor¬ 
responded to roughly lOMEdd accretion rate. All the simulations re¬ 
sulted in optically thick and geometrically thick turbulent accretion 
disks. In the course of this study we reached a number of conclu¬ 
sions: 

(i) Photosphere: Only for accretion rates M < lOMEdd does the 
photosphere extend down to the horizon, i.e., only for such low ac¬ 
cretion rates will an observer at infinity viewing down the axis be 
able to see the innermost region of the accretion disk. For larger ac¬ 
cretion rates, an on-axis observer will only observe a photosphere 
which is located relatively far from the black hole. In fact, even 
for relatively low accretion rates M < lOMEdd^ the optically thin 
region is limited to polar angles 6 < 15°. For larger inclinations, 
an observer would only see the photosphere of the accretion disk 
located at large distance from the BH, often outside the computa¬ 
tional box of the simulations. The limited box size and duration of 
the simulations do not allow us to study the radiation coming from 
such distant photospheres in any detail. We expect, however, that 
the spectrum of any such radiation will be relatively soft, because 
of the large distance from the BH, and that the observed isotropic- 
equivalent luminosity will not significantly exceed the Eddington 
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value - a super-Eddington flux in an optically thick wind causes 
acceleration of the wind and reduces the radiative flux towards the 
Eddington limit. 

(ii) Stagnation radius: Eor the simulations with the lowest ac¬ 
cretion rates (M < lOMEdd). the stagnation radius in the funnel, tq, 
which separates an inner region where gas falls into the BH from an 
outer region where gas flows out, is located near the BH {tq < 10). 
In the case of simulation rOll (with a spinning BH), where the 
extraction of the rotational energy of the BH increases the energy 
flux through the funnel, the stagnation radius is very close to the 
horizon. Eor the run with the largest accretion rate (r003), the stag¬ 
nation radius moves signiflcantly out - all gas in the funnel within 
ro ^ 50 falls on the BH. This is a result of the increased optical 
depth of the gas, which traps photons and effectively suppresses 
the outflow of radiation up to this radius. 

(iii) Total luminosity: All the simulations with non-rotating BHs 
show the same total efflciency of roughly 3%Mc^, approximately 
a factor of two less than the efflciency of a standard thin accretion 
disk. We obtain similar ratio of efficiencies for the simulation with 
a rotating BH (a^ = 0.7), for which the measured total efflciency 
was ~ ^%Mc^. These efficiencies were calculated from the total 
luminosity in all forms of energy: radiative, kinetic, magnetic, ther¬ 
mal and binding (only rest mass energy was not included). The total 
luminosity is thus fundamental and represents the total energy de¬ 
posited at “infinity” (the interstellar medium). There is no unique 
way of decomposing the above total luminosity into its constituent 
parts because of the limited size of the inflow/outflow equilibrium 
region in the simulations. When measured at the limit of inflow- 
outflow equilibrium (r ^ 25 in the disk interior), energy balance is 
dominated by the binding energy flux. This energy must be trans¬ 
formed into other forms of energy before reaching infinity (where 
the binding energy is zero). 

(iv) Radiative luminosity: Radiative luminosity can be measured 
reliably only inside the optically thin funnel region near the axis; 
here, radiation is decoupled from gas, and whatever radiation is 
flowing outward is guaranteed to reach the observer. However, only 
the simulations with the lowest accretion rates show an optically 
thin region inside the computational domain at all. Even in these 
cases, the optically thin radiative luminosity increases with radius 
(Tablebecause of radiation flowing into the funnel from the disk 
at larger radii. To obtain a useful estimate of the net radiative lumi¬ 
nosity from the funnel, we would have to simulate accretion flows 
in much bigger computational box and for a much longer time. 
This is presently impractical. The radiative luminosities measured 
at radius r - 250 in the optically thin and outflowing regions are 
quite low. Eor accretion rates near lOMEdd and a non-rotating BH, 
only ~ 20% of the total liberated energy of 3%Mc^ (mentioned 
above) comes out as optically thin radiation escaping through the 
funnel. The kinetic luminosity exceeds the radiative luminosity for 
the largest accretion rate considered, where the coupling between 
radiation and gas in the funnel is strongest, resulting in radiative 
acceleration of gas ( [S^dowski & Narayan||2015b| ). The luminosi¬ 
ties are signiflcantly larger for the simulation with a rotating BH 
because the BH spin provides an extra source of energy. 

(v) Beaming: We have confirmed that radiation is beamed along 
the polar axis and the radiative fluxes here can be signiflcantly 
super-Eddington. Eor the fiducial accretion rate of lOMEdd. the sim¬ 
ulations with non-rotating and spinning BHs show radiative fluxes 
of 20 and lOOFEdd^ respectively, on the axis. In contrast, the kinetic 
energy is not beamed on the axis but peaks (Eig.[^ either in a coni¬ 
cal shell (for lower accretion rates) or in the wind region (for larger 
accretion rates). 


(vi) Photon trapping: We compared the total flux of radiation 
with the diffusive flux and showed that, for the accretion rates con¬ 
sidered, advection of radiation (photon trapping) dominates over 
diffusive transport deep inside the disk. Eor M ^ lOMEdd. the ad- 
vective flux is more than 30 times stronger than the diffusive flux 
near the equatorial plane, and this ratio increases by at least an order 
of magnitude when the accretion rate grows to ~ 176MEdd. which 
is explained by the increased optical depth. Closer to the disk sur¬ 
face or funnel boundary, the diffusive flux of radiation dominates, 
which ultimately contributes to the optically thin radiation escap¬ 
ing through the funnel. This flux also helps to drive gas away from 
the disk surface into the funnel. 

(vii) Radiation transport An analysis of our simulations shows 
that radiation transport is dominated by advective and diffusive 
transport. We do not see a significant component of turbulent ra¬ 
diative transport. 

(viii) Advection factor: Eor the fiducial model (M = lOMEdd) we 
estimated the fraction of photons generated in the disk that manage 
to escape vertically and enter the optically thin funnel region. At 
radius r = 10, only ~ 20% of photons leave the disk while ~ 80% 
end up in the BH. Corresponding numbers at radius r - 25 are 
approximately 45% and 55%, respectively. This suggests that the 
effective trapping radius, where half the photons escape from the 
disk and half fall into the BH, is located at r > 30 for lOMEdd- 

(ix) Variability: We attempted to extract from the simulations 
(frequency integrated) lightcurves as seen by observers at infinity. 
Because of the limited size of the equilibrium region, and the fact 
that the photosphere is located close to the axis even for the low¬ 
est accretion rate considered, we were able to obtain robust light 
curves only for radiation escaping along the axis, and that too only 
for the fiducial model. Extracting variability information from sim¬ 
ulations of optically thick and geometrically very thick disks is, 
and will remain to be, challenging. One may expect that thin disks, 
with photospheres much closer to the equatorial plane, will be more 
amenable to study. 

(x) Impact of the BH mass: We performed one simulation 
(r6)20) with BH mass lOOOM©, which had the same Eddington- 
scaled mass accretion rate as the fiducial model (rOOl). The prop¬ 
erties of the two simulations, after scaling down the latter by appro¬ 
priate factors, were quite similar; in fact, it was difficult to distin¬ 
guish the two (compare Eig.|^. This is because electron scattering 
opacity dominates in both simulations, and under these conditions 
the accretion equations scale very simply with BH mass. Eor much 
larger BH masses (corresponding to the SMBH regime), the ab¬ 
sorption opacity will no longer be negligible and this will break 
exact scaling. 

(xi) Axisymmetric simulations: We compared the properties 
of the fiducial three-dimensional simulation (rOOl) with a two- 
dimensional axisymmetric simulation (d300) that made use of an 
artificial magnetic dynamo and accreted at roughly the same rate. 
We showed that the time-averaged properties of the two simula¬ 
tions were remarkably similar. Because of the signiflcantly lower 
computational cost, we were able to run the axisymmetric simu¬ 
lation for a much longer time and obtained a signiflcantly larger 
equilibrium region (extending up to r ~ 100 instead of only r ~ 25 
in the case of the 3D model). Thus, two-dimensional axisymmet¬ 
ric simulations (with magnetic dynamo) are a cheap and promising 
method for running models for long times and thereby extending 
the range of radii over which one obtains useful information. How¬ 
ever, this must be done with caution. While time-averged quanti¬ 
ties appear to be reliable, we note that the temporal properties of 
the 2D and 3D runs were different. In particular, because of the 
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lack of non-axisymmetric modes in the axisymmetric simulation, 
it showed much larger variability in the radiative flux flowing out 
along the axis. 


8.1 Comparison with other studies 


Our study is the first to explore the parameter space relevant to 
radiation-dominated BH accretion disks using three-dimensional, 
general relativistic, radiation-MHD simulations. The properties of 
the simulated super-critical disks described in this work are quali¬ 
tatively in agreement with previously published global simulations. 
In particular, significant photon trapping was identified already 
by|Ohsug a & Min eshige| ( |2007|), and confir med more recently by 
[Sqdowski et al. ( |2014| ) and McKinney et al.| ( [20T4) . The total lumi¬ 
nosities given in the present work are also close to those obtained 
in the latter studies. The same is true for radiation beaming and the 
radiative luminosity of the funnel. 

There are a number of differences between our results and 
those described in |Jiang et al.| ( |2014b] ). The model presented in the 
latter work, which had an accretion rate of ~ 15MEdd. can be di¬ 
rectly compared to our fiducial model rOOl. pTang et al.| ( |2014b| ) 
found a very different spatial distribution of gas compared to our 
run; their disk is strongly concentrated at the equatorial plane 
whereas our disk is not. Perhaps because of this, they obtained 
a powerful radiative flux from the innermost (r < 10) region; in 
fact, nearly all of their radiative luminosity, which is of the order 
of 5 %Mc 2, comes from such small radii. Our simulations show a 
similar total efficiency (3%Mc^), but only a small fraction of the 
energy escapes as radiation at small radii. 

|Jiang et ^ ( |2014b| ) used a Newtonian code with cylindrical 
coordinates (and a “cylindrical BH” with radius r = A). Their simu¬ 
lation was intialized with a different and more strongly magnetized 
torus compared to ours. They implemented a sophisticated radia¬ 
tion transfer solver which evolved in real time a number of specific 
intensities. In comparison, our simulations were done with the Ml 
closure scheme, which means that we evolved only four indepen¬ 
dent radiation quantities in each cell. On the other hand, | Jiang et al.| 
( |2014b| ), for the sake of performance, had to make some approxi¬ 
mations when treating the interactions between gas and radiatioij^ 
Which of these factors is responsible for the large discrepancy be¬ 
tween their study and ours is presently unclear. 
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^ These approximations were avoided by |Ohsuga & TakahashiH2015] , who 
implemented sub-stepping in their implicit solver, but otherwise used the 
same ideas as in | Jiang et al.| ( |2014^ . 
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